
> ################################################################
> # Replication File: "Backscratching in Banks: Political Cycles in Bank Manager Appointments."
> # Author: Jonas Markgraf
> #
> # Main Analysis contains (in chronological order as appearing in paper):
> ## -- replication code for all empirical analyses and plots in the paper
> #
> # R Version: 3.6.0 (2019-04-26)
> # platform:  x86_64-apple-darwin15.6.0
> ################################################################
> 
> rm(list = ls())

> # load packages
> library(eha)

> library(survival)

> library(flexsurv)

> library(texreg)  

> library(dplyr)

> library(simPH)

> # load data
> 
> ## data for Spanish savings banks
> cajas.boards <- read.table(file = "cajas_dta.txt", header = T)  

> ###############################################################################
> ############## MAIN RESULTS ###################################################
> ###############################################################################
> 
> # **Table 1: Main Results**:  -----------------------------------------------
> 
> ## Column 1
> fit.main1 <- coxph(Surv(bank.yrsInOffice, bank.event) ~  region.yrsInGov,
+   method = "efron",
+   data = cajas.boards, subset = bank.chair == 1)

> summary(fit.main1)
Call:
coxph(formula = Surv(bank.yrsInOffice, bank.event) ~ region.yrsInGov, 
    data = cajas.boards, subset = bank.chair == 1, method = "efron")

  n= 1259, number of events= 234 

                    coef exp(coef) se(coef)      z Pr(>|z|)    
region.yrsInGov -0.05595   0.94559  0.01127 -4.963 6.94e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
region.yrsInGov    0.9456      1.058    0.9249    0.9667

Concordance= 0.604  (se = 0.019 )
Likelihood ratio test= 28.07  on 1 df,   p=1e-07
Wald test            = 24.63  on 1 df,   p=7e-07
Score (logrank) test = 25.39  on 1 df,   p=5e-07


> ## Column 2
> fit.main2 <- coxph(Surv(bank.yrsInOffice, bank.event) ~  region.postElec,
+   method = "efron",
+   data = cajas.boards, subset = bank.chair == 1)

> summary(fit.main2)
Call:
coxph(formula = Surv(bank.yrsInOffice, bank.event) ~ region.postElec, 
    data = cajas.boards, subset = bank.chair == 1, method = "efron")

  n= 1259, number of events= 234 

                  coef exp(coef) se(coef)     z Pr(>|z|)   
region.postElec 0.3494    1.4182   0.1338 2.611  0.00902 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
region.postElec     1.418     0.7051     1.091     1.843

Concordance= 0.549  (se = 0.018 )
Likelihood ratio test= 6.95  on 1 df,   p=0.008
Wald test            = 6.82  on 1 df,   p=0.009
Score (logrank) test = 6.89  on 1 df,   p=0.009


> ##Column 3
> fit.main3 <- coxph(Surv(bank.yrsInOffice, bank.event) ~ region.yrsInGov +
+     region.postElec +
+     region.coal + 
+     bank.polVote +
+     bank.roaLag1 +
+     as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice +
+     frailty(bank.id, distribution = "gaussian"),
+   method = "efron",
+   data = cajas.boards, subset = bank.chair == 1)

> summary(fit.main3)
Call:
coxph(formula = Surv(bank.yrsInOffice, bank.event) ~ region.yrsInGov + 
    region.postElec + region.coal + bank.polVote + bank.roaLag1 + 
    as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice + 
    frailty(bank.id, distribution = "gaussian"), data = cajas.boards, 
    subset = bank.chair == 1, method = "efron")

  n= 1259, number of events= 234 

                          coef     se(coef) se2     Chisq DF    p      
region.yrsInGov           -0.05776 0.01575  0.01487 13.44  1.00 2.5e-04
region.postElec            0.28108 0.13964  0.13908  4.05  1.00 4.4e-02
region.coal               -0.20405 0.17897  0.17299  1.30  1.00 2.5e-01
bank.polVote              -0.03099 0.90536  0.75864  0.00  1.00 9.7e-01
bank.roaLag1               0.03954 0.14522  0.13819  0.07  1.00 7.9e-01
as.factor(region.govParty -0.08520 0.31997  0.28459  0.07  1.00 7.9e-01
as.factor(region.govParty  0.80198 0.28840  0.26518  7.73  1.00 5.4e-03
frailty(bank.id, distribu                           67.07 36.43 1.5e-03
as.factor(region.govParty -7.18812 0.86614  0.86595 68.87  1.00 1.0e-16
as.factor(region.govParty -7.21841 0.86577  0.86558 69.51  1.00 7.6e-17
as.factor(region.govParty -7.30733 0.86602  0.86585 71.20  1.00 3.2e-17

                          exp(coef) exp(-coef) lower .95 upper .95
region.yrsInGov           0.9438785     1.0595 0.9151807  0.973476
region.postElec           1.3245562     0.7550 1.0074191  1.741529
region.coal               0.8154216     1.2264 0.5741750  1.158031
bank.polVote              0.9694863     1.0315 0.1643974  5.717267
bank.roaLag1              1.0403301     0.9612 0.7826391  1.382868
as.factor(region.govParty 0.9183316     1.0889 0.4905013  1.719329
as.factor(region.govParty 2.2299466     0.4484 1.2670902  3.924474
as.factor(region.govParty 0.0007555  1323.6158 0.0001383  0.004126
as.factor(region.govParty 0.0007330  1364.3203 0.0001343  0.004000
as.factor(region.govParty 0.0006706  1491.1833 0.0001228  0.003661

Iterations: 6 outer, 120 Newton-Raphson
     Variance of random effect= 0.3520437 
Degrees of freedom for terms=  0.9  1.0  0.9  0.7  0.9  1.6 36.4  2.8 
Concordance= 0.946  (se = 0.946 )
Likelihood ratio test= 814.3  on 45.3 df,   p=<2e-16


> ## Column 4
> fit.main4a <- coxph(Surv(bank.yrsInOffice, bank.event) ~ region.postElec +
+     region.coal:region.postElec +
+     region.coal + 
+     bank.polVote +
+     bank.roaLag1 +
+     as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice +
+     frailty(bank.id, distribution = "gaussian"),
+   method = "efron",
+   data = cajas.boards, subset = bank.chair == 1 & region.govChange == 1)

> summary(fit.main4a)
Call:
coxph(formula = Surv(bank.yrsInOffice, bank.event) ~ region.postElec + 
    region.coal:region.postElec + region.coal + bank.polVote + 
    bank.roaLag1 + as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice + 
    frailty(bank.id, distribution = "gaussian"), data = cajas.boards, 
    subset = bank.chair == 1 & region.govChange == 1, method = "efron")

  n= 445, number of events= 94 

                          coef    se(coef) se2    Chisq DF p      
region.postElec            1.5228 0.4637   0.4622 10.78  1 1.0e-03
region.coal                1.0937 0.4620   0.4574  5.61  1 1.8e-02
bank.polVote               1.0261 1.2142   1.1034  0.71  1 4.0e-01
bank.roaLag1              -0.2720 0.2443   0.2369  1.24  1 2.7e-01
as.factor(region.govParty  0.6194 0.4530   0.4322  1.87  1 1.7e-01
as.factor(region.govParty  0.1791 0.5670   0.5496  0.10  1 7.5e-01
frailty(bank.id, distribu                         13.20 10 2.1e-01
region.postElec:region.co -1.7452 0.5357   0.5336 10.61  1 1.1e-03
as.factor(region.govParty -8.2219 1.0658   1.0655 59.51  1 1.2e-14
as.factor(region.govParty -8.4778 1.0640   1.0636 63.49  1 1.6e-15
as.factor(region.govParty -8.5291 1.0646   1.0642 64.18  1 1.1e-15

                          exp(coef) exp(-coef) lower .95 upper .95
region.postElec           4.5850192     0.2181 1.848e+00 11.377400
region.coal               2.9854025     0.3350 1.207e+00  7.382883
bank.polVote              2.7902019     0.3584 2.583e-01 30.140539
bank.roaLag1              0.7618337     1.3126 4.720e-01  1.229770
as.factor(region.govParty 1.8578533     0.5383 7.646e-01  4.514532
as.factor(region.govParty 1.1960866     0.8361 3.936e-01  3.634343
region.postElec:region.co 0.1746033     5.7273 6.110e-02  0.498947
as.factor(region.govParty 0.0002687  3721.4671 3.327e-05  0.002170
as.factor(region.govParty 0.0002080  4807.0598 2.585e-05  0.001674
as.factor(region.govParty 0.0001976  5060.1299 2.453e-05  0.001592

Iterations: 6 outer, 120 Newton-Raphson
     Variance of random effect= 0.1682414 
Degrees of freedom for terms=  1.0  1.0  0.8  0.9  1.8 10.0  1.0  2.8 
Concordance= 0.93  (se = 0.93 )
Likelihood ratio test= 306.4  on 19.37 df,   p=<2e-16


> ## Column 5
> fit.main4b <- coxph(Surv(bank.yrsInOffice, bank.event) ~ region.postElec +
+     region.coal:region.postElec +
+     region.coal + 
+     bank.polVote +
+     bank.roaLag1 +
+     as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice +
+     frailty(bank.id, distribution = "gaussian"),
+   method = "efron",
+   data = cajas.boards, subset = bank.chair == 1 & region.govChange == 0)

> summary(fit.main4b)
Call:
coxph(formula = Surv(bank.yrsInOffice, bank.event) ~ region.postElec + 
    region.coal:region.postElec + region.coal + bank.polVote + 
    bank.roaLag1 + as.factor(region.govParty) + as.factor(region.govParty):bank.yrsInOffice + 
    frailty(bank.id, distribution = "gaussian"), data = cajas.boards, 
    subset = bank.chair == 1 & region.govChange == 0, method = "efron")

  n= 814, number of events= 140 

                          coef     se(coef) se2    Chisq DF    p      
region.postElec            0.36048 0.1920   0.1910  3.53  1.00 6.0e-02
region.coal               -0.34419 0.4488   0.4376  0.59  1.00 4.4e-01
bank.polVote              -1.42344 1.2115   1.0824  1.38  1.00 2.4e-01
bank.roaLag1               0.08905 0.1888   0.1793  0.22  1.00 6.4e-01
as.factor(region.govParty -0.63922 0.5074   0.4495  1.59  1.00 2.1e-01
as.factor(region.govParty  1.14472 0.4139   0.3686  7.65  1.00 5.7e-03
frailty(bank.id, distribu                          33.16 21.87 5.8e-02
region.postElec:region.co  0.59112 0.5531   0.5516  1.14  1.00 2.9e-01
as.factor(region.govParty -6.46281 1.4924   1.4923 18.75  1.00 1.5e-05
as.factor(region.govParty -6.43082 1.4921   1.4920 18.57  1.00 1.6e-05
as.factor(region.govParty -6.55022 1.4920   1.4920 19.27  1.00 1.1e-05

                          exp(coef) exp(-coef) lower .95 upper .95
region.postElec            1.434018     0.6973 9.844e-01   2.08907
region.coal                0.708793     1.4108 2.941e-01   1.70814
bank.polVote               0.240884     4.1514 2.242e-02   2.58862
bank.roaLag1               1.093133     0.9148 7.550e-01   1.58263
as.factor(region.govParty  0.527703     1.8950 1.952e-01   1.42667
as.factor(region.govParty  3.141552     0.3183 1.396e+00   7.07011
region.postElec:region.co  1.806008     0.5537 6.109e-01   5.33934
as.factor(region.govParty  0.001560   640.8576 8.373e-05   0.02908
as.factor(region.govParty  0.001611   620.6828 8.650e-05   0.03001
as.factor(region.govParty  0.001430   699.3977 7.678e-05   0.02663

Iterations: 7 outer, 140 Newton-Raphson
     Variance of random effect= 0.2783625 
Degrees of freedom for terms=  1.0  1.0  0.8  0.9  1.6 21.9  1.0  2.9 
Concordance= 0.955  (se = 0.955 )
Likelihood ratio test= 504.9  on 30.98 df,   p=<2e-16


> ## print all results
> screenreg(list(fit.main1, fit.main2, fit.main3, fit.main4a, fit.main4b),
+   caption = "Cox PH Model", caption.above = T,
+   label="T:findingsChair",
+   custom.coef.names = c("Years in Government"
+     , "Post Election"
+     , "Coalition"
+     , "Return on Assets$_{t-1}$"
+     , "Public Sector Vote Share"
+     , "Party: Other", "Party: Socialist"
+     , "Post Election * Coalition"),
+   custom.model.names = c("Raw Model 1", "Raw Model 2", "Full Model", "After Gov't Change", "No Gov't Change"),
+   include.rsquared=F, include.maxrs=F, include.missings=F, include.zph=F,
+   reorder.coef = c(1:2,8,3:7),
+   stars = c(.01,.05,.1),
+   omit.coef = c("bank.yrsInOffice"))

=====================================================================================================
                           Raw Model 1  Raw Model 2  Full Model   After Gov't Change  No Gov't Change
-----------------------------------------------------------------------------------------------------
Years in Government          -0.06 ***                 -0.06 ***                                     
                             (0.01)                    (0.02)                                        
Post Election                              0.35 ***     0.28 **     1.52 ***             0.36 *      
                                          (0.13)       (0.14)      (0.46)               (0.19)       
Post Election * Coalition                                          -1.75 ***             0.59        
                                                                   (0.54)               (0.55)       
Coalition                                              -0.20        1.09 **             -0.34        
                                                       (0.18)      (0.46)               (0.45)       
Return on Assets$_{t-1}$                               -0.03        1.03                -1.42        
                                                       (0.91)      (1.21)               (1.21)       
Public Sector Vote Share                                0.04       -0.27                 0.09        
                                                       (0.15)      (0.24)               (0.19)       
Party: Other                                           -0.09        0.62                -0.64        
                                                       (0.32)      (0.45)               (0.51)       
Party: Socialist                                        0.80 ***    0.18                 1.14 ***    
                                                       (0.29)      (0.57)               (0.41)       
-----------------------------------------------------------------------------------------------------
AIC                        2928.33      2949.45      2230.69      705.32              1225.16        
Num. events                 234          234          234          94                  140           
Num. obs.                  1259         1259         1259         445                  814           
=====================================================================================================
*** p < 0.01, ** p < 0.05, * p < 0.1

> ## **Figure 1**: Plot effect of time in government ----------------------------
> 
> ## simulate values from Table 1, model 3, using 'simPH' package
> plot.region.yrsInGov <- coxsimLinear(fit.main3
+   , b = "region.yrsInGov"
+   , Xj = seq(1, 28, by = 1)) # for year 1--28
All Xl set to 0.

> ## plot
> fig1 <- simGG(plot.region.yrsInGov
+   , xlab = "Years in Government"
+   , lcolour = "black"
+   , pcolour = "grey"
+   , alpha = .5
+   , type = "points")

> pdf("figures/effect_govChangeYr.pdf")

> fig1
`geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

> dev.off()
RStudioGD 
        2 

> # **Figure 2**: Plot interactions from Table 1 -------------------------------
> 
> ## prepare coefficients from Table 1, model 4, for plotting
> 
> ### extract values
> fit.main4a$coefficients
                                        region.postElec                                             region.coal 
                                              1.5227943                                               1.0937346 
                                           bank.polVote                                            bank.roaLag1 
                                              1.0261140                                              -0.2720269 
                       as.factor(region.govParty)Others                     as.factor(region.govParty)Socialist 
                                              0.6194217                                               0.1790551 
                            region.postElec:region.coal as.factor(region.govParty)Conservative:bank.yrsInOffice 
                                             -1.7452388                                              -8.2218732 
      as.factor(region.govParty)Others:bank.yrsInOffice    as.factor(region.govParty)Socialist:bank.yrsInOffice 
                                             -8.4778409                                              -8.5291474 

> relevant.coefs1 <- as.data.frame(coef(fit.main4a)[c(1:2,7)])

> ### obtain full interaction
> relevant.coefs1[3,] <- sum(coef(fit.main4a)[c(1:2, 7)])  

> ### calculate standard error for interaction
> se_interact1 <- sqrt(fit.main4a$var2[1,1] + fit.main4a$var2[2,2] + 2*fit.main4a$var2[1,2])

> ### extract values for 90% confidence interval
> relevant.cis1 <- as.data.frame(confint(fit.main4a, level = .9)[c(1:2,7), c(1:2)]) 

> ### estimate 90% confidence interval for full interaction
> relevant.cis1[3,1] <- relevant.coefs1[3,] - 1.645*se_interact1

> relevant.cis1[3,2] <- relevant.coefs1[3,] + 1.645*se_interact1

> ### bind everything together and clean
> exp.int1 <- cbind(relevant.coefs1, relevant.cis1)

> colnames(exp.int1)[1] <- "est"

> ### drop main effect for coalition
> exp.int1 <- exp.int1[-2,]

> exp.int1
                                  est        5 %     95 %
region.postElec             1.5227943  0.7600762 2.285512
region.postElec:region.coal 0.8712901 -0.5654770 2.308057

> ## prepare coefficients from Table 1, model 5, for plotting
> 
> fit.main4b$coefficients
                                        region.postElec                                             region.coal 
                                             0.36048029                                             -0.34419137 
                                           bank.polVote                                            bank.roaLag1 
                                            -1.42343783                                              0.08904789 
                       as.factor(region.govParty)Others                     as.factor(region.govParty)Socialist 
                                            -0.63922130                                              1.14471687 
                            region.postElec:region.coal as.factor(region.govParty)Conservative:bank.yrsInOffice 
                                             0.59111887                                             -6.46280726 
      as.factor(region.govParty)Others:bank.yrsInOffice    as.factor(region.govParty)Socialist:bank.yrsInOffice 
                                            -6.43082021                                             -6.55021951 

> relevant.coefs2 <- as.data.frame(coef(fit.main4b)[c(1:2,7)])

> relevant.coefs2[3,] <- sum(coef(fit.main4b)[c(1:2, 7)])

> se_interact2 <- sqrt(fit.main4b$var2[1,1] + fit.main4b$var2[2,2] + 2*fit.main4b$var2[1,2])

> relevant.cis2 <- as.data.frame(confint(fit.main4b, level = .9)[c(1:2,7), c(1:2)])

> relevant.cis2[3,1] <- relevant.coefs2[3,] - 1.645*se_interact2

> relevant.cis2[3,2] <- relevant.coefs2[3,] + 1.645*se_interact2

> exp.int2 <- cbind(relevant.coefs2, relevant.cis2)

> colnames(exp.int2)[1] <- "est"

> exp.int2 <- exp.int2[-2,]

> exp.int2
                                  est         5 %     95 %
region.postElec             0.3604803  0.04473162 0.676229
region.postElec:region.coal 0.6074078 -0.23913694 1.453953

> ## plot the obtained estimates for Table 1, model 4 and 5
> 
> pdf("figures/interaction_log-hazardratio.pdf")

> par (mar=c(4,6,0,0))

> plot (c(1,3), c(floor(min(exp.int1[2])), ceiling(max(exp.int1[3]))), type="n", bty="n"
+   , ylab=""
+   , xlab=""
+   , axes=F)

> y_axis <- seq(floor(min(exp.int1[2])), ceiling(max(exp.int1[3])), by=1)

> axis (2, at = y_axis, labels = parse(text = paste("e^", y_axis)), las = 1)

> mtext (text="Relative Hazard", side=2, line=2.5, cex = 2)

> mtext (text=c("Single Party", "Coalition"), side=1, at=c(1.5,2.5), line=0.5, cex = 2)

> abline(h=0, lty =2)

> points (xy.coords(c(1.4,2.4), exp.int1[,1])
+   , pch=19, col=c("black","black"), cex=2)

> segments (x0=c(1.4,2.4), x1=c(1.4,2.4), y0=exp.int1[,2], y1=exp.int1[,3]
+   , lwd=3, col=c("black", "black"))

> points (xy.coords(c(1.6,2.6), exp.int2[,1])
+   , pch=19, col=c("grey","grey"), cex=2)

> segments (x0=c(1.6,2.6), x1=c(1.6,2.6), y0=exp.int2[,2], y1=exp.int2[,3]
+   , lwd=3, col=c("grey", "grey"))

> legend("topright", legend = c("New Government", "Reelected Government"), col = c("black", "grey"), lty = 1
+   , pch = c(19, 19), box.lty=0, cex = 1.5)

> dev.off()
RStudioGD 
        2 
